module ObsData

use Index
use strutils
use parameter
use SysUtils 
use HTypes

implicit none


contains


function GetDataSourceTechName(DataSourceType, DataSourceIndex)
   implicit none
   integer, intent(in) :: DataSourceType, DataSourceIndex
   character(len = cTechNameLength) :: GetDataSourceTechName
   select case(DataSourceType)
      case(odCanum)
         GetDataSourceTechName = 'CANUM'
      case(odWeCa)
         GetDataSourceTechName = 'WECA'
      case(odWeSt)
         GetDataSourceTechName = 'WEST'
      case(odMatProp)
         GetDataSourceTechName = 'MATPROP'
      case(odFleet)
         GetDataSourceTechName = 'FLEET ' // lrtrim(inttostr(DataSourceIndex + 1))
      case(odSSB)
         GetDataSourceTechName = 'SSB ' // lrtrim(inttostr(DataSourceIndex + 1))
      case(odBiomass)
         GetDataSourceTechName = 'BMASS ' // lrtrim(inttostr(DataSourceIndex + 1))
      case default
         GetDataSourceTechName = ''
   end select
end function GetDataSourceTechName 
   
   
function GetObsData(Opt, ObservedDataType, ObservedDataIndex)
   implicit none
   type(TOptimisation), intent(in) :: Opt
	integer, intent(in) :: ObservedDataType, ObservedDataIndex
	type(TObservedData), pointer :: GetObsData 
	GetObsData => Opt%ObservedData(ObservedDataType)%Elm(ObservedDataIndex)
end function GetObsData


function GetNumObservedDataElms(Opt, ObservedDataType)
   implicit none
   type(TOptimisation), intent(in) :: Opt
   integer, intent(in) :: ObservedDataType
   integer :: GetNumObservedDataElms
   GetNumObservedDataElms = 0
   if (associated(Opt%ObservedData(ObservedDataType)%Elm)) then
      GetNumObservedDataElms = size(Opt%ObservedData(ObservedDataType)%Elm)
   endif
end function GetNumObservedDataElms


function GetObsDataValue(Opt, ObservedDataType, ObservedDataIndex, iYear, iAge)
! Purpose: Returns the observed data value from population survey,  year and age index.
! In:      ObsDataType -  the observed data type.
!          iSurveyIndex - Survey index into population model.
!          iYear - Year index into population model.
!          iAge - Age index into population model.
! Return   -1 if not defined, otherwise return the observed value
! Note:    iPopYear and iPopAge are converted into iObsYear and iObsAge.
!          If not found in observed data, return -1 which means 'not used'
	implicit none
	double precision :: GetObsDataValue
   type(TOptimisation), intent(in) :: Opt
	integer, intent(in) :: ObservedDataType
	integer, intent(in) :: ObservedDataIndex, iYear, iAge
	integer :: iObsYear, iObsAge 
	type(TObservedData), pointer :: pVar 
	GetObsDataValue = -1
	
	pVar => GetObsData(Opt, ObservedDataType, ObservedDataIndex)
   !if (.not. associated(pVar)) return
   iObsYear = iYear + Opt%IndexInput%FirstYear - pVar%FirstYear
   if (pVar%FirstAge == -1) then
      iObsAge = 0
   else
	   iObsAge = iAge + Opt%IndexInput%FirstAge - pVar%FirstAge
   endif
	if ((0 <= iObsYear) .and. (iObsYear <= pVar%NumYears - 1) .and. &
		(0 <= iObsAge) .and. (iObsAge <= pVar%NumAges - 1)) then
		GetObsDataValue = pVar%Obs(iObsYear, iObsAge)
	end if
end function GetObsDataValue

 
function GammaCategoryFromDataSource(DataSourceType)
   implicit none
   integer :: DataSourceType
   integer :: GammaCategoryFromDataSource
   select case(DataSourceType)
      case(odCanum)
         GammaCategoryFromDataSource = cGamma_Canum
      case(odFleet)
         GammaCategoryFromDataSource = cGamma_Fleet
      case(odSSB)
         GammaCategoryFromDataSource = cGamma_SSB
      case(odBiomass)
         GammaCategoryFromDataSource = cGamma_BM
   end select
end function GammaCategoryFromDataSource

 
function GetFRef(Opt, iYear)
   implicit none
   double precision :: GetFRef
   type(TOptimisation), intent(in) :: Opt
   double precision :: NRef
   integer :: iFirstRef, iLastRef, iAge, iYear
   iFirstRef = Opt%IndexInput%FirstRefAge - Opt%IndexInput%FirstAge
   iLastRef = Opt%IndexInput%LastRefAge - Opt%IndexInput%FirstAge
   GetFRef = 0
   NRef = 0
   do iAge = iFirstRef, iLastRef
      if (Opt%PopModel%Fay(iYear, iAge) > -1) then
         GetFRef = GetFRef + Opt%PopModel%Fay(iYear, iAge)
         NRef = NRef + 1
      endif
   end do
   if (0 < NRef) then
      GetFRef = GetFRef / NRef
   end if
end function GetFRef


subroutine SetResDataValue(Opt, ObservedDataType, ObservedDataIndex, iYear, iAge, ModVal, PhiVal, PhiWeight, ResVal)
! Purpose: Set residual value from population survey,  year and age index.
! In:      ObsDataType -  the observed data type.
!          iSurveyIndex - Survey index into population model.
!          iYear - Year index into population model.
!          iAge - Age index into population model.
! Return   -1 if not defined, otherwise return the observed value
! Note:    iPopYear and iPopAge are converted into iObsYear and iObsAge.
!          If not found in observed data, return -1 which means 'not used'
	implicit none
   type(TOptimisation), intent(in) :: Opt
	integer, intent(in) :: ObservedDataType
	integer, intent(in) :: ObservedDataIndex, iYear, iAge
	double precision, intent(in) :: ModVal, PhiVal, PhiWeight, ResVal
	double precision :: CblFacVal, CblPowVal, F, M, Z, C, FProp, MProp, NFromDS, ObsVal, FyRef
	integer :: iObsYear, iObsAge 
	type(TObservedData), pointer :: pVar 
	type(TDataSource), pointer :: pDataSource
	
	pVar => GetObsData(Opt, ObservedDataType, ObservedDataIndex)
   !if (.not. associated(pVar)) return
   iObsYear = iYear + Opt%IndexInput%FirstYear - pVar%FirstYear
	iObsAge = iAge + Opt%IndexInput%FirstAge - pVar%FirstAge
	if ((0 <= iObsYear) .and. (iObsYear <= pVar%NumYears - 1) .and. &
		(0 <= iObsAge) .and. (iObsAge <= pVar%NumAges - 1)) then
		pVar%ModVal(iObsYear, iObsAge) = ModVal
		pVar%ResU(iObsYear, iObsAge) = ResVal
		if(PhiWeight /= 0.0d0) then
		   pVar%ResArray(pVar%NumPositiveUnweightedResiduals) = ResVal
		   pVar%NumPositiveUnweightedResiduals = pVar%NumPositiveUnweightedResiduals + 1
		endif
		pVar%ResW(iObsYear, iObsAge) = ResVal * max(0.0d0, PhiWeight)
		pVar%Phi(iObsYear, iObsAge) = PhiVal
		pVar%PhiWeight(iObsYear, iObsAge) = max(0.0d0, PhiWeight)
	   M = PValue(Opt, cNatMort, 0, iYear, iAge) 
	   F = Opt%PopModel%Fay(iYear, iAge)
	   Z = M + F
	   select case(ObservedDataType)
	      case (odCanum)
	         !CblVal = PValue(Opt, cSelection, 0, iYear, iAge) 
	         !Use normalized selection instead.
	         FyRef = GetFRef(Opt, iYear)
	         if ((F > 0.0d0) .and. (FyRef > 0.0d0)) then
	            CblFacVal = F / FyRef
	         endif
	         CblPowVal = 1.0d0
	         ! Calculate N from catch
	      case (odFleet)
	         CblFacVal = PValue(Opt, cCBL_Fac_Fleet, ObservedDataIndex, iYear, iAge) 
	         CblPowVal = PValue(Opt, cCBL_Pow_Fleet, ObservedDataIndex, iYear, iAge) 
	      case (odSSB)
	         CblFacVal = PValue(Opt, cCBL_Fac_SSB, ObservedDataIndex, iYear, 0) 
	         CblPowVal = PValue(Opt, cCBL_Pow_SSB, ObservedDataIndex, iYear, 0) 
	      case (odBiomass)
	         CblFacVal = PValue(Opt, cCBL_Fac_BM, ObservedDataIndex, iYear, 0) 
	         CblPowVal = PValue(Opt, cCBL_Pow_BM, ObservedDataIndex, iYear, 0) 
	      case default
	         CblFacVal = 0.0d0
	         CblPowVal = 0.0d0
	   end select
   	NFromDS = 0.0d0
	   select case(ObservedDataType)
	      case (odCanum)
            C = GetObsDataValue(Opt, odCaNum, 0, iYear, iAge)
            if(F > 0.0d0) then
   	         NFromDS = C * Z / F / (1 - exp(-Z))
   	      endif
   	   case (odFleet)
            pDataSource => GetDataSource(Opt, ObservedDataType, ObservedDataIndex)
            FProp = 0
            MProp = 0
            if (associated(pDataSource)) then
               FProp = pDataSource%FProp
               MProp = pDataSource%MProp
            endif
            ObsVal = GetObsDataValue(Opt, ObservedDataType, ObservedDataIndex, iYear, iAge)
            if ((CblFacVal > 0.0d0) .and. (CblPowVal > 0.0d0) .and. (F >= 0.0d0) .and. (ObsVal >= 0.0d0)) then
               NFromDS = ((ObsVal / CblFacVal)**(1 / CblPowVal)) * exp(M * MProp + F * FProp)
            endif
	      case default
	         NFromDS = -1
	   end select
		pVar%CblFac(iObsYear, iObsAge) = CblFacVal
		pVar%CblPow(iObsYear, iObsAge) = CblPowVal
		pVar%NFromDS(iObsYear, iObsAge) = NFromDS
	end if
end subroutine SetResDataValue


subroutine CreateNewObsData(Opt)
   implicit none
   type(TOptimisation), intent(inout) :: Opt
   integer :: v, N, ioerr, I
   type(TObservedData), pointer :: pVar

   ! Reallocate observed data when it is read from file
   call DeallocateObsData(Opt)
   allocate(Opt%ObservedData(odFirst : odLast))
   do v = odFirst, odLast
      N = 1
      select case(v)
         case (odFleet, odFleetWeight)
            N = NumFleetIndices(Opt%IndexInput)
         case (odSSB, odSSBWeight)
            N = NumSSBIndices(Opt%IndexInput)
         case (odBiomass, odBiomassWeight)
            N = NumBiomassIndices(Opt%IndexInput)
      end select 
      if (N == 0) cycle
      allocate(Opt%ObservedData(v)%Elm(0 : N - 1))
      do I = 0, N - 1
         pVar => Opt%ObservedData(v)%Elm(I)
         pVar%Title = "Observed data"
         pVar%Variable = -1
         pVar%FirstYear = 2000
         pVar%NumYears = 6
         pVar%FirstAge = 1
         select case(v)
            case(odCaton, odSSB, odSSBWeight, odBiomass, odBiomassWeight)
               pVar%FirstAge = -1
               pVar%NumAges = 1
            case default
               pVar%NumAges = 10
         end select
         call AllocateObservedData(pVar)
         select case(v)
            case(odCanumWeight, odFleetWeight, odSSBWeight, odBiomassWeight)
               pVar%Obs(:, :) = 1.0d0
            case default
               pVar%Obs(:, :) = 0.0d0
         end select
         call GetNumberFormatInMatrix(pVar%Obs, pVar%NumDig, pVar%NumDec)
      enddo
   end do
end subroutine CreateNewObsData


function ReadObsData(Opt)
   implicit none
   type(TOptimisation), intent(inout) :: Opt
   integer :: v, N, I
   type(TObservedData), pointer :: pVar
   logical :: ReadObsData
   ReadObsData = .false.

   ! Reallocate observed data when it is read from file
   call DeallocateObsData(Opt)
   allocate(Opt%ObservedData(odFirst : odLast))

   do v = odFirst, odLast
      N = 1
      select case(v)
         case (odFleet, odFleetWeight)
            N = NumFleetIndices(Opt%IndexInput)
         case (odSSB, odSSBWeight)
            N = NumSSBIndices(Opt%IndexInput)
         case (odBiomass, odBiomassWeight)
            N = NumBiomassIndices(Opt%IndexInput)
      end select 
      if (0 < N) then
         allocate(Opt%ObservedData(v)%Elm(0 : N - 1))
         select case(v)
            case (odCanumWeight, odFleetWeight, odSSBWeight, odBiomassWeight)
               if(.not. OpenFile(Opt, cFileUnit, GetWorkFileName(Opt, &
                  lrtrim(Opt%IndexInput%DataFiles(v))), .false.)) return
            case default
               if(.not. OpenFile(Opt, cFileUnit, GetDataFileName(Opt, &
                  lrtrim(Opt%IndexInput%DataFiles(v))), .false.)) return
         end select
         
         select case(v)
            case (odFleet, odFleetWeight)
               if(.not. ReadAgeStructuredIndices(Opt, v)) return
            case (odSSB, odSSBWeight, odBiomass, odBiomassWeight)
               if(.not. ReadNonAgeStructuredIndices(Opt, v)) return
            case default
               if(.not. ReadObservedData(Opt, v)) return
         end select 
         close(cFileUnit)
         do I = 0, size(Opt%ObservedData(v)%Elm) - 1
            pVar => Opt%ObservedData(v)%Elm(I)
            call GetNumberFormatInMatrix(pVar%Obs, pVar%NumDig, pVar%NumDec)
         enddo
      endif
   end do
   ReadObsData = .true.
end function ReadObsData


subroutine DeAllocateObservedData(pVar)
   implicit none 
   type(TObservedData), pointer :: pVar

   if (associated(pVar)) then
      if (associated(pVar%Obs)) then
         deallocate(pVar%Obs)
      end if
      if (associated(pVar%ModVal)) then
         deallocate(pVar%ModVal)
      end if
      if (associated(pVar%ResU)) then
         deallocate(pVar%ResU)
      end if
      if (associated(pVar%ResW)) then
         deallocate(pVar%ResW)
      end if
      if (associated(pVar%Phi)) then
         deallocate(pVar%Phi)
      end if
      if (associated(pVar%PhiWeight)) then
         deallocate(pVar%PhiWeight)
      end if
      if (associated(pVar%PhiArray)) then
         deallocate(pVar%PhiArray)
      end if
      if (associated(pVar%CblFac)) then
         deallocate(pVar%CblFac)
      end if
      if (associated(pVar%CblPow)) then
         deallocate(pVar%CblPow)
      end if
      if (associated(pVar%NFromDS)) then
         deallocate(pVar%NFromDS)
      end if
   endif
end subroutine DeAllocateObservedData


subroutine AllocateObservedData(pVar)
   implicit none 
   type(TObservedData), pointer :: pVar
   ! Ensure it is deallocated
   call DeAllocateObservedData(pVar)
   
   allocate(pVar%Obs(0 : (pVar%NumYears - 1), 0 : (pVar%NumAges - 1)))
   allocate(pVar%ModVal(0 : (pVar%NumYears - 1), 0 : (pVar%NumAges - 1)))
   allocate(pVar%ResU(0 : (pVar%NumYears - 1), 0 : (pVar%NumAges - 1)))
   allocate(pVar%ResW(0 : (pVar%NumYears - 1), 0 : (pVar%NumAges - 1)))
   allocate(pVar%Phi(0 : (pVar%NumYears - 1), 0 : (pVar%NumAges - 1)))
   allocate(pVar%PhiWeight(0 : (pVar%NumYears - 1), 0 : (pVar%NumAges - 1)))
   allocate(pVar%PhiArray(0 : (pVar%NumYears * pVar%NumAges - 1)))
   allocate(pVar%ResArray(0 : (pVar%NumYears * pVar%NumAges - 1)))
   allocate(pVar%CblFac(0 : (pVar%NumYears - 1), 0 : (pVar%NumAges - 1)))
   allocate(pVar%CblPow(0 : (pVar%NumYears - 1), 0 : (pVar%NumAges - 1)))
   allocate(pVar%NFromDS(0 : (pVar%NumYears - 1), 0 : (pVar%NumAges - 1)))
   ! Reset values
   pVar%Obs(0 : pVar%NumYears - 1, 0 : pVar%NumAges - 1) = 0.0d0
   pVar%ModVal(0 : pVar%NumYears - 1, 0 : pVar%NumAges - 1) = 0.0d0
   pVar%ResU(0 : pVar%NumYears - 1, 0 : pVar%NumAges - 1) = 0.0d0
   pVar%ResW(0 : pVar%NumYears - 1, 0 : pVar%NumAges - 1) = 0.0d0
   pVar%Phi(0 : pVar%NumYears - 1, 0 : pVar%NumAges - 1) = 0.0d0
   pVar%PhiWeight(0 : pVar%NumYears - 1, 0 : pVar%NumAges - 1) = 0.0d0
   pVar%PhiArray(0 : pVar%NumYears * pVar%NumAges - 1) = 0.0d0
   pVar%ResArray(0 : pVar%NumYears * pVar%NumAges - 1) = 0.0d0
   pVar%CblFac(0 : pVar%NumYears - 1, 0 : pVar%NumAges - 1) = 0.0d0
   pVar%CblPow(0 : pVar%NumYears - 1, 0 : pVar%NumAges - 1) = 0.0d0
   pVar%NFromDS(0 : pVar%NumYears - 1, 0 : pVar%NumAges - 1) = 0.0d0
   pVar%NumPositiveUnweightedResiduals = 0
end subroutine AllocateObservedData

 
function ReadObservedData(Opt, v)
   implicit none
   type(TOptimisation), intent(in) :: Opt
   integer, intent(in) :: v
   integer :: iDummy, LastYear, LastAge, A, Y
   type(TObservedData), pointer :: pVar
   character(len=30) :: ElmName
   logical :: ReadObservedData
   ReadObservedData = .false.
   
   pVar => Opt%ObservedData(v)%Elm(0)

   ! Read the stock title

   if(.true.) then
      if(.not. ReadStrFromFile(Opt, cFileUnit, pVar%Title)) return
   	if(.not. ReadIFromFile(Opt, cFileUnit, iDummy, -1, -1, "Not used")) return
   	if(.not. ReadIFromFile(Opt, cFileUnit, pVar%Variable, -1, -1, "Variable type", .true.)) return
   	if(.not. ReadIFromFile(Opt, cFileUnit, pVar%FirstYear, 0, -1, "First year")) return
   	if(.not. ReadIFromFile(Opt, cFileUnit, LastYear, 0, -1, "Last year", .true.)) return
   	if(.not. CheckRange(Opt, pVar%FirstYear, LastYear, "First year", "Last year")) return
   	if(.not. ReadIFromFile(Opt, cFileUnit, pVar%FirstAge, -1, -1, "First age")) return
   	if(.not. ReadIFromFile(Opt, cFileUnit, LastAge, -1, -1, "Last age", .true.)) return
   	if(.not. CheckRange(Opt, pVar%FirstAge, LastAge, "First age", "Last age")) return
   	if(.not. ReadIFromFile(Opt, cFileUnit, iDummy, -1, -1, "Not used", .true.)) return
   else
      read(cFileUnit, '(A)') pVar%Title
      read(cFileUnit, *) iDummy, pVar%Variable
      read(cFileUnit, *) pVar%FirstYear, LastYear, pVar%FirstAge, LastAge, iDummy
   endif
   pVar%NumYears = LastYear - pVar%FirstYear + 1
   pVar%NumAges = LastAge - pVar%FirstAge + 1
   if (v == odCaTon) then
      pVar%FirstAge = -1
      pVar%NumAges = 1
   endif
   call AllocateObservedData(pVar)
   do Y = 0, pVar%NumYears - 1
      if(.true.) then
         do A = 0, pVar%NumAges - 1
            ElmName = "Year:" // lrtrim(inttostr(pVar%FirstYear + Y)) // ", Age:" // &
      	      lrtrim(inttostr(pVar%FirstAge + A))
      	   if(.not. ReadFFromFile(Opt, cFileUnit, pVar%Obs(Y, A), -1.0d0, -1.0d0, &
      	      ElmName, pVar%NumAges - 1 == A)) return
         enddo
      else
         read(cFileUnit, *) pVar%Obs(Y, :)
      endif
   enddo
   ReadObservedData = .true.
end function ReadObservedData


function ReadAgeStructuredIndices(Opt, v)
   implicit none
   type(TOptimisation), intent(inout) :: Opt
   integer :: I, Y, A, NumOfFleetsPlus100, LastAge, LastYear
   double precision :: effort, CPUE, TimeBegin, TimeEnd
   integer :: sex
   integer, intent(in) :: v
   type(TObservedData), pointer :: pVar
   character(len=80) :: FileTitle
   character(len=30) :: ElmName
   logical :: ReadAgeStructuredIndices
   ReadAgeStructuredIndices = .false.
   
   if(.true.) then
      if(.not. ReadStrFromFile(Opt, cFileUnit, FileTitle)) return
   	if(.not. ReadIFromFile(Opt, cFileUnit, NumOfFleetsPlus100, 100, -1, "Number of fleets + 100", .true.)) return
      do I = 0, size(Opt%ObservedData(v)%Elm) - 1
         pVar => Opt%ObservedData(v)%Elm(I)
         if(.not. ReadStrFromFile(Opt, cFileUnit, pVar%Title)) then
            ! EOF - allocate empty fleet
		      pVar%NumYears = 1  
		      pVar%NumAges = 1  
		      pVar%FirstYear = 2000
		      pVar%FirstAge = 1
            call AllocateObservedData(pVar)
            Opt%ErrorStr = ""
            cycle
         endif
         if(.not. ReadIFromFile(Opt, cFileUnit, pVar%FirstYear, 0, -1, "First year")) return
         if(.not. ReadIFromFile(Opt, cFileUnit, LastYear, 0, -1, "Last year", .true.)) return
      	if(.not. CheckRange(Opt, pVar%FirstYear, LastYear, "First year", "Last year")) return
         if(.not. ReadIFromFile(Opt, cFileUnit, sex, -1, -1, "Sex")) return
         if(.not. ReadFFromFile(Opt, cFileUnit, effort, -1.0d0, -1.0d0, "Effort")) return
         if(.not. ReadFFromFile(Opt, cFileUnit, TimeBegin, -1.0d0, -1.0d0, "TimeBegin")) return
         if(.not. ReadFFromFile(Opt, cFileUnit, TimeEnd, -1.0d0, -1.0d0, "TimeEnd", .true.)) return
         if(.not. ReadIFromFile(Opt, cFileUnit, pVar%FirstAge, 0, -1, "First age")) return
         if(.not. ReadIFromFile(Opt, cFileUnit, LastAge, 0, -1, "Last age", .true.)) return
		   pVar%NumYears = LastYear - pVar%FirstYear + 1  
		   pVar%NumAges = LastAge - pVar%FirstAge + 1  
         call AllocateObservedData(pVar)
         do Y = 0, pVar%NumYears - 1
            if(.true.) then
               if(.not. ReadFFromFile(Opt, cFileUnit, CPUE, -1.0d0, -1.0d0, "CPUE")) return
               do A = 0, pVar%NumAges - 1
                  ElmName = "Year:" // lrtrim(inttostr(pVar%FirstYear + Y)) // ", Age:" // &
      	            lrtrim(inttostr(pVar%FirstAge + A))
      	         if(.not. ReadFFromFile(Opt, cFileUnit, pVar%Obs(Y, A), -1.0d0, -1.0d0, &
      	            ElmName, pVar%NumAges - 1 == A)) return
               enddo
            else
               read(cFileUnit, *) pVar%Obs(Y, :)
            endif
         enddo
      end do 
   else
      read(cFileUnit, '(A)') FileTitle
      read(cFileUnit, *) NumOfFleetsPlus100
      do I = 0, size(Opt%ObservedData(v)%Elm) - 1
         pVar => Opt%ObservedData(v)%Elm(I)
         read(cFileUnit, '(A)') pVar%Title
         read(cFileUnit, *) pVar%FirstYear, LastYear
         read(cFileUnit, *) sex, effort, TimeBegin, TimeEnd
         read(cFileUnit, *) pVar%FirstAge, LastAge
		   pVar%NumYears = LastYear - pVar%FirstYear + 1  
		   pVar%NumAges = LastAge - pVar%FirstAge + 1  
         call AllocateObservedData(pVar)
         do Y = 0, pVar%NumYears - 1
            read(cFileUnit, *) CPUE, pVar%Obs(Y, :)
         enddo
      end do 
   endif

   ReadAgeStructuredIndices = .true.
end function ReadAgeStructuredIndices


function ReadNonAgeStructuredIndices(Opt, v)
   implicit none
   type(TOptimisation), intent(in) :: Opt
   integer :: NumInd, Y, Yr, I
   double precision :: VPA
   integer, intent(in) :: v
   type(TObservedData), pointer :: pVar
   double precision, pointer :: SourceObs(:, :)
   character(len=40) :: FileTitle
   character(len = 8) :: YrLbl, VPALbl
   character(len = 8), pointer :: Labels(:)
   integer :: ObsData, NumIndInFile, NumYearsInFile, VPALoc
   logical :: ReadNonAgeStructuredIndices
   ReadNonAgeStructuredIndices = .false.
   
	NumInd = size(Opt%ObservedData(v)%Elm)
	read(cFileUnit, '(A)') FileTitle
	read(cFileUnit, *) NumIndInFile, NumYearsInFile, VPALoc
	allocate(Labels(0 : NumInd - 1))
	allocate(SourceObs(0 : NumYearsInFile - 1, 0 : NumIndInFile - 1));
	read(cFileUnit, *) YrLbl, VPALbl, Labels
	do Y = 0, NumYearsInFile - 1 
		read(cFileUnit, *) Yr, VPA, SourceObs(Y, :)
		do I = 0, NumInd - 1
      	pVar => Opt%ObservedData(v)%Elm(I)
			if (Y == 0) then
				pVar%Title = Labels(I)
				pVar%FirstYear = Yr
				pVar%NumYears = NumYearsInFile
				pVar%FirstAge = 0
				pVar%NumAges = 1
            call AllocateObservedData(pVar)
			end if
			pVar%Obs(Y, 0) = SourceObs(Y, I)
		end do
	end do
	deallocate(Labels)
	deallocate(SourceObs)
	ReadNonAgeStructuredIndices = .true.
end function ReadNonAgeStructuredIndices


subroutine WriteObsData(Opt)
   implicit none
   type(TOptimisation), intent(inout) :: Opt
   integer :: v, N, ioerr
   do v = odFirst, odLast
      N = 1
      select case(v)
         case (odFleet, odFleetWeight)
            N = NumFleetIndices(Opt%IndexInput)
         case (odSSB, odSSBWeight)
            N = NumSSBIndices(Opt%IndexInput)
         case (odBiomass, odBiomassWeight)
            N = NumBiomassIndices(Opt%IndexInput)
      end select 
      if (0 < N) then
         select case(v)
            case (odCanumWeight, odFleetWeight, odSSBWeight, odBiomassWeight)
               open(cFileUnit, file = GetWorkFileName(Opt, lrtrim(Opt%IndexInput%DataFiles(v))), &
                  action = 'write', iostat = ioerr, status = 'replace')
            case default
               open(cFileUnit, file = GetDataFileName(Opt, lrtrim(Opt%IndexInput%DataFiles(v))), &
                  action = 'write', iostat = ioerr, status = 'replace')
         end select
         if (io_no_error .ne. ioerr) then
            call MessageOut(Opt, "Error writing :" // lrtrim(Opt%IndexInput%DataFiles(v)))
         endif
         select case(v)
            case (odFleet, odFleetWeight)
               call WriteAgeStructuredIndices(Opt, v)
            case (odSSB, odSSBWeight, odBiomass, odBiomassWeight)
               call WriteNonAgeStructuredIndices(Opt, v)
            case default
               call WriteObservedData(Opt, v)
         end select 
         close(cFileUnit)
      endif
   end do
end subroutine WriteObsData


subroutine WriteObservedData(Opt, v)
   implicit none
   type(TOptimisation), intent(in) :: Opt
   integer, intent(in) :: v
   integer :: Misc, DimType, LastYear, LastAge, Y, A
   type(TObservedData), pointer :: pVar
   pVar => Opt%ObservedData(v)%Elm(0)
   call GetNumberFormatInMatrix(pVar%Obs, pVar%NumDig, pVar%NumDec)

   ! Write the stock title
	LastYear = pVar%NumYears + pVar%FirstYear - 1  
	LastAge = pVar%FirstAge + pVar%NumAges - 1    
   if (v == odCaTon) then
      LastAge = -1
   endif
   ! Bug corrected. Write Misc=1(allways) and DimType=1(array)
   ! Note DimType 2 for  vector maybe. Test with R.
   Misc = 1
   DimType = 1
   write(cFileUnit, '(A)') lrtrim(pVar%Title)
   write(cFileUnit, *) Misc, pVar%Variable
   write(cFileUnit, *) pVar%FirstYear, LastYear
   write(cFileUnit, *) pVar%FirstAge, LastAge
   write(cFileUnit, *) DimType
   do Y = 0, pVar%NumYears - 1
      write(cFileUnit, "(1000" // lrtrim(GetOutFormat(pVar%NumDig, pVar%NumDec)) // ")") pVar%Obs(Y, :)
   enddo
end subroutine WriteObservedData


function GetOutFormatColWidth(NumDig, NumDec)
   implicit none
   integer :: GetOutFormatColWidth, NumDig, NumDec
   GetOutFormatColWidth = 3 + NumDig + 1 + NumDec
end function GetOutFormatColWidth


function GetOutFormat(NumDig, NumDec)
   implicit none
   character(len = 15) :: GetOutFormat
   integer :: NumDig, NumDec
   GetOutFormat = "F" // lrtrim(inttostr(GetOutFormatColWidth(NumDig, NumDec))) // "." &
      // lrtrim(inttostr(NumDec)) 
end function GetOutFormat


subroutine WriteAgeStructuredIndices(Opt, v)
   implicit none
   type(TOptimisation), intent(in) :: Opt
   integer :: I, Y, A, NumOfFleetsPlus100, LastAge, LastYear
   integer, intent(in) :: v
   type(TObservedData), pointer :: pVar
   write(cFileUnit, '(A)') "Age structured tuning indices"
   NumOfFleetsPlus100 = size(Opt%ObservedData(v)%Elm) + 100 ! Lowestoft standard (?).
   write(cFileUnit, *) NumOfFleetsPlus100

   do I = 0, size(Opt%ObservedData(v)%Elm) - 1
      pVar => Opt%ObservedData(v)%Elm(I)
      call GetNumberFormatInMatrix(pVar%Obs, pVar%NumDig, pVar%NumDec)
      write(cFileUnit, '(A)') lrtrim(pVar%Title)
		LastYear = pVar%NumYears + pVar%FirstYear - 1  
      write(cFileUnit, *) pVar%FirstYear, LastYear
      write(cFileUnit, *) -1, -1, -1, -1  
		LastAge = pVar%FirstAge + pVar%NumAges - 1    
      write(cFileUnit, *) pVar%FirstAge, LastAge
      do Y = 0, pVar%NumYears - 1
         write(cFileUnit, "(I, 1000" // lrtrim(GetOutFormat(pVar%NumDig, pVar%NumDec)) // ")") -1, pVar%Obs(Y, :)
      enddo
   end do 
end subroutine WriteAgeStructuredIndices


subroutine WriteNonAgeStructuredIndices(Opt, v)
   implicit none
   type(TOptimisation), intent(in) :: Opt
   integer :: NumInd, Y, Yr, I
   integer :: VPA
   integer, intent(in) :: v
   type(TObservedData), pointer :: pVar
   double precision, pointer :: SourceObs(:, :)
   character(len=40) :: FileTitle
   character(len = 8), pointer :: Labels(:)
   integer :: ObsData, NumYearsInFile, VPALoc
	NumInd = size(Opt%ObservedData(v)%Elm)
	FileTitle = "SSB indices"
	write(cFileUnit, '(A)') FileTitle
	NumYearsInFile = Opt%ObservedData(v)%Elm(0)%NumYears
	VPALoc = 2
	VPA = -1
	write(cFileUnit, *) NumInd, NumYearsInFile, VPALoc
	allocate(Labels(0 : NumInd - 1))
	allocate(SourceObs(0 : NumYearsInFile - 1, 0 : NumInd - 1));
	do I = 0, NumInd - 1
	   Labels(I) = '   INDEX' // lrtrim(inttostr(I))
	enddo
	write(cFileUnit, *) "YEAR", "   VPA", Labels
	do Y = 0, NumYearsInFile - 1 
		do I = 0, NumInd - 1
      	pVar => Opt%ObservedData(v)%Elm(I)
         call GetNumberFormatInMatrix(pVar%Obs, pVar%NumDig, pVar%NumDec)
   	   Yr = pVar%FirstYear + Y
   	   SourceObs(Y, I) =  pVar%Obs(Y, 0)
		end do
      write(cFileUnit, "(I4,2x, I2, 2x, 1000" // lrtrim(GetOutFormat(pVar%NumDig, pVar%NumDec)) // ")") Yr, VPA, SourceObs(Y, :)
	end do
	deallocate(Labels)
	deallocate(SourceObs)
end subroutine WriteNonAgeStructuredIndices


subroutine DeallocateObsData(Opt)
!   Purpose: This procedure frees allocated memory
!   History: When        Who        What
!            02-06-07    Skalevik   Created
   implicit none
   type(TOptimisation), intent(inout) :: Opt
   type(TObservedData), pointer :: pVar
   integer :: v, I, Num
   if (associated(Opt%ObservedData)) then
      do v = odFirst, odLast
         if (associated(Opt%ObservedData(v)%Elm)) then
            Num = size(Opt%ObservedData(v)%Elm)
            if (0 < Num) then
               do I = 0, Num - 1
                  pVar => Opt%ObservedData(v)%Elm(I)
                  call DeAllocateObservedData(pVar)
               end do
               deallocate(Opt%ObservedData(v)%Elm)
            end if
         end if
      end do
      deallocate(Opt%ObservedData)
   end if
end subroutine DeallocateObsData


end module ObsData
